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C"» ' ABSTRACT 

We extend to large scales a method proposed in previous work that reconstructs non- 
parametrically the primordial power spectrum from cosmic microwave background 
data at high resolution. The improvement is necessary to account for the non- 
gaussianity of the Wilkinson Microwave Anisotropy Probe (WMAP) likelihood due 
primarily to cosmic variance. We assume the concordance ACDM cosmology, utilise a 
smoothing prior and perform Monte Carlo simulations around an initial power spec- 
. trum that is scale-free and with spectral index n s = 0.97, very close to the concor- 

' dance spectrum. The horizon scale for the model we are considering corresponds to 

l/" - ) | the wavenumber kh = 4.52 x 10 _4 Mpc . We find some evidence for the presence 

of features and we quantify the probabilities of exceeding the observed deviations in 
WMAP data with respect to the fiducial models. We detect the following marginal 
departures from a scale-free (spectral index n s = 0.97) initial spectrum: a cut-off at 
0.0001 < k < 0.001 Mpc" 1 at 79.5% (92%), a dip at 0.001 < k < 0.003 Mpc" 1 at 
87.2% (98%) and a bump at 0.003 < k < 0.004 Mpc" 1 at 90.3% (55.5%) confidence 
+3 , level. 

These frequentist confidence levels are calculated by integrating over the distribu- 
tion of the Monte Carlo reconstructions built around the fiducial models. The frequen- 
tist analysis finds the low k cutoff of the estimated power spectrum to be about 2.5c 
away from the n s = 0.97 model, while in the Bayesian analysis the model is about 
1.5<7 away from the estimated spectrum. (The a's are different for the two different 
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1 INTRODUCTION 

The spectacular results provided by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) experiment (Bennett 
et al. 2003) offer the possibility of testing the nature of the 
initial conditions that seeded the structure of the Universe. 
The inflationary paradigm puts forward a convincing pic- 
ture that explains the formation of the first seeds. However 
inflation has not been yet incorporated successfully in the 
framework of fundamental physics. Therefore hints from ob- 
servations may prove useful in discerning such an important 
connection. 

An increasing number of studies have been devoted 
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to the initial power spectrum recon struction, due to it s 
potential im p ortanc e. For example in lSper gel et alJ 12003): 
iPeiris et al.l (120031) : iBridle et all i2003l) : iBarger et alJ 
J2003h : Mukheriee & Wang (2003); Mukherjee & Wang 
(2005) a shape parametrized a priori is used. Model 
independent methods that use so-called non para- 
metric algorithms ha ve also been proposed e .g. in 



■ Gawiser fc Silkl <ll99cl) ; iTegniark fc ZaldarriaiL 
Matsumiva ct al j2002ll2003h:lKogo et alJ j2004aUbl 
Shafieloo fc Souradeed i2004T ) : iHu fc Okamotol 



2002) 
2005) 



2004) 



iHannestadl <2004lh iTocchini-Valentini et all (120051): Leach 
(2005) ; Sealfon et al. (2005). In ITocchini-Valentini et alJ 
(2005) we proposed a method based on regularized least 
squares to non-parametrically reconstruct under a general 
smoothing assumption the primordial power spectrum 
from CMB temperature data in the wave number range 
0.01 < k < O.lftMpc" 1 , with h standing for the Hubble 
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1 Mpc 1 . Here we extend 
proposing the necessary 



constant in units of 100 Km sec - 
the analysis to large scales by 
modifications dictated by the non-gaussianity of the data 
errors caused by cosmic variance. While we fix the late-time 
cosmological parameters to the best fit values obtained from 
the combined WMAP and Sloan Digital Sky Survey (SDSS) 
data (Tegmark et al. 2004), our method has the virtue of 
allowing an estimate of the power spectrum, at the highest 
resolution to date, that permits a refined local search for 
deviations from a scale free power spectrum. The horizon 
scale for the model we are considering translates into the 
wavenumber kh ~ 2-k/tio = 4.52 x 10 -4 Mpc -1 , where 770 is 
the conformal time at present. A high resolution reconstruc- 
tion at large scales (using the Rich ardson-Lucy algorithm) 
has b een presented previously by IShafieloo fc Souradeed 
(2004), that used binned data after the hexadecapole. 
We propose here a different deconvolution method that 
starts from the unbinned WMAP data and provides for the 
first time an estimate of the errors in the reconstruction. 
Deconvolving WMAP data we find indications for the 
following three features that deviate from a fiducial scale 
free (spectral index n 3 = 0.97) initial spectrum: a cut-off 
at 0.0001 < k < 0.001 Mpc -1 at 79.5% (92%), a dip at 
0.001 < k < 0.003 Mpc -1 at 87.2% (98%) and a bump 
at 0.003 < k < 0.004 Mpc -1 at 90.3% (55.5%). Similar 
featur es were also detected in IShafieloo fc Souradeed 
i2004Tl an d suggestions of a si milar c utoff were a l so pro - 
vided bv iMukheriee fc Wand (I2003D: iBridle et al.l J2003>) : 
IContaldi et alJ J2003I) : ICline et alJ J2003I) . The cutoff effect 
is due mainly to the low large scale power at the quadrupole 
and to a minor extent to the octopole present in WMAP 
data (Bennett et al. 2003) and, since it is intriguingly 
present also in the COBE DMR data (Bennett et al. 1996) , 
certainly deserves some consideration. 

We wish to emphasize that our approach is very similar 
in spirit to traditional image reconstruction methods and 
that our aim is not to search for the minimal parametric 
description of the power spectrum. Given a dense enough 
discretization to closely approximate the integral in Eq.Q 
below, the data automatically select the degree of smooth- 
ing necessary to allow the highest resolution in the recon- 
structed signal compatible with the actual noise level, being 
the result of a tradeoff between resolution and noise filter- 
ing. Strictly speaking we are using a large number of band 
amplitudes as parameters, however since they correspond 
to a dense grid in wavenumber space that translates in a 
good approximation of the continuum, we call our method 
non-parametric. As will become apparent later, our method 
is capable of estimating simultaneously the large number of 
amplitude parameters. Widely used multi-parameter esti- 
mation methods such as Monte Carlo Markov Chain would, 
in their simplest implementations, find some difficulties in 
evaluating hundreds of parameters. After having tested the 
reliability of the method through Monte Carlo simulations, 
the reconstruction can be compared with simple parametric 
shapes to look for local deviations, that we estimated with 
two types of errors, as detailed below. 

In Section II we introduce in detail the method; in Sec- 
tion III we show the outcome of some tests on mock data 
that allow us to proceed with confidence to the reconstruc- 
tion from the actual WMAP temperature data, our main 



result. We conclude with a d iscussion of our results in Sec- 
tion IV. 



2 ESTIMATION OF THE POWER SPECTRUM 
2.1 Introduction to the Method 

The problem at hand is how to estimate the primordial 
power spectrum P(k) from a CMB observational data base 
which is compressed into the form of the angular power spec- 
trum, Ci. 

We start first with a pedagogical example where the 
cosmic variance is neglected, a situation applicable in the 
large / case. 

The temperature CMB angular power spectrum is re- 
lated to the primordial power spectrum, P(k), by a convo- 
lution with a window function that depends on the cosmo- 
logical parameters: 

C e = 4ir J A 2 e (k)P°(k)^~Y^WuXi, (1) 

where the last term represents a numerical approximation 
to the integral obtained through a fine discretization and 
Xi = P(ki)/ki. We fix the units in the transfer function 
in order to express the primordial fluctuation spectrum x, 
correspondi ng to the comoving curvature perturbation spec- 
trum A|, in I Verde etaD <2003ft . in units of 2.95 x 10" 9 . 

Let us suppose that our scope is to estimate a solution 
for the linear problem 



d = Wx + n, 



(2) 



where d is the CMB vector data, n is the noise vector, char- 
acterized by zero mean and covariance matrix 



N 



(3) 



and x is the vector of the unknowns. The CMB measure- 
ments are given by the angular power spectrum evaluated at 
the multipoles £, W is the transfer function calculated nu- 
merically with a Boltzmann code such as CMBFAST (Sel- 
jak & Zaldarriaga 1996), and x represents the initial power 
spectrum, as defined above. Cosmic variance is neglected 
here and N is independent of x. 

Consider a Gaussian likelihood, 



C tx exp 



-i(Wa; - d)'l\l -1 (WiE - d) 







= exp 





, (4) 



where we are ignoring an overall factor that does not depend 
on the model and Eq.(|lJ is implicitly defining S". Our pur- 
pose is to find an estimate for the vector x, given the data, 
their covariance matrix and the window function. 
The solution to the minimization problem 



d 



Si WHO, 



(•>) 



maximizes the likelihood and provides a minimum variance 
estimator that is unbiased. The solution is found to be 



x = [W*N -1 W] 1 W*N -1 d, 

and its covariance matrix is given by 

e = rw f N _1 wi _1 . 



(6) 



(7) 
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However in this context it is not possible to compute the in- 
verse matrix defined above, essentially because the transfer 
function W does not represent a one-to-one relation connect- 
ing the data space to the solution space, since the projection 
on the last scattering surface and smearing effects (including 
gravitational lensing, that is however not considered in the 
present work) destroy part of the information (Tegmark & 
Zaldarriaga 2002). Note that the formal solution provided 
by Eqs.JIjJ and J7J is valid only under the assumption that 
N is independent of x. 

Clearly to solve the problem, some external informa- 
tion must be injected in the form of a prior. In order to do 
this in a meaningful way, the problem must be rephrased in 
Bayesian terms. We are interested in evaluating the proba- 
bility P(x\d) of having a primordial power spectrum given 
the data. Bayes theorem states that: 



P(x\d) = 



P{d\x)P{x) 
P(d) 



(8) 



where P(d\x) is the likelihood of the data given a model 
for the power spectrum, P(x) is the prior assumed on the 
power spectrum and P(d) is the probability of having the 
data. In what follows we will fix the late-ti me cosmological 
param eters to the best fit values found by iTeemark et alJ 
(2004) when they assumed a featurless power law for the 
initial power spectrum. 

We will consider a prior that favours smooth solutions 
and a useful choice is to adopt a prior of the form 



P(x) oc exp ( — ex t L t Lx) = exp (— eR) . 



(9) 



with L representing an approximation to the first deriva- 
tive operator and e is a constant parameter. Eq.© defines 
R as the regularizing operator whose strength depends on 
e. In other words, the logarithm of the prior, namely the 
regularization function, approximates a penalty term that 
corresponds to the squared norm of the first derivative of 
the solution. After neglecting the probability of the data, 
that does not depend on the initial power spectrum, the 
minimization problem generalizes to 



{S'[x] + eR[x]} =0. 



(10) 



d_ 

Ox 

The parameter e weighs the amount of smoothing that the 
solution is requested to have. In other words, it is possible to 
find a solution, provided it is smoothed to a certain degree. 
Under the assumption that the covariance matrix N does 
not depend on the theoretical power spectrum, the solution 
can be written as 

x = [W*N : W + elA] _1 W'N^d, (11) 

with the now computable covariance matrix 

E = [W'N 1 ^/ + eL'L] _1 . (12) 

The dependence of the solution on the smoothness param- 
eter e appears clearly in the previous equations. It is useful 
to view some asymptotic limits of the solution. When the 
parameter e tends to zero we are back to the first case of no 
prior and the solution is going to be impossible to find for 
the reasons exposed above. If e is taken very large and the 
smoothing term dominates over the first term, the solution 
will be in the form that minimizes the first derivative, that 
is, a constant vector. 



Clearly, a compromise is needed when the parameter e 
has to be chosen. A reasonable choice is to select an e such 
that 



S'[x] = N d 



(13) 



with Nd representing the number of data points. In essence, 
this was the method adopted in our first reconstruction pa- 
per, iTocchini-Valentini et alJ (120051) . where, for high enough 
multipoles, cosmic variance could be considered negligible 
with respect to experimental noise and a Gaussian shape 
for the likelihood of the data given a model that was a safe 
choice, as will also be confirmed below. 

Clearly the algorithm can now be cast as a classical 
Lagrange multiplier problem. Namely the equations to be 
solved are 



S'[x] = N d 
and 

^{S'[x] + eR[x}}=0. 



(14) 



(15) 



2.2 Extension to Large Scales 

Since we are interested in an estimate of the power spectrum 
at large scales, corresponding to large angle CMB data, the 
simple Gaussian likelihood of the previous section can no 
longer be used. The first non-Gaussian modifications to get 
a reliable fit for the CMB likelihood and the need for them 
were proposed in Bond, Jaffe, Knox (1998 and 2000). For the 
WMAP experiment, the likelihood is approximated in detail 
by the product of a lognormal and a Gaussian distributions 
(Verde et al. 2003): 

- 21n£ = ^{z th - z d ) t Q- 1 (z th - z d ) + 



1 



Wx-dflVT^Waj-d), 



(16) 



with z th = ln(Wa; + Af) and z d = \n(d + Af). Here 
Q-) = (Wa; + Af) e N~} (\Nx + N) t , . 



(17) 

The noise vector A/" was computed by the WMAP team 
(Verde et al. (2003)) through calibration with Monte Carlo 
simulations. It is important to note that the matrix N de- 
pends on the angular power spectrum Wa;, due to cosmic 
variance and in addition the galaxy cut provokes some cor- 
relation due to mode mixing. 

Since the primordial power spectrum reconstruction on 
large scales is very sensitive to the quadrupole likelihood 
function, it is relevant to be sure that foreground contamina- 
tion has been properly taken into account. There have been 
various post- WMAP papers in which alternative methods 
to clean foregrounds have been presented (e.g. in Tegmark 
et al. (2003), Efstathiou (2004) and Slosar et al. (2004)). 
We will discuss later the possible effects of alternative fore- 
ground removal methods. 

Even if in this work we chose to use the Verde et 
al.(2003) likelihhood, our method can be generalised pro- 
vided an analytical fit to the likelihood is given. 

Application of Bayes theorem and the requirement that 
the estimated power spectrum should correspond to the 
maximum of the a posteriori probability reduces the mini- 
mization problem to 
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A {5 + ei?} = , 



where 

s 



z d ) + 



2/ th d\tg-\ — l/ th 

-(z - z ) Q (z 



r (Wx - d)'N _1 (Wj; - d). 



under the constraint of 



S[x] = N d . 



(18) 



(19) 



(20) 



Here the covariance matrix N has an explicit dependence on 
x to account for the cosmic variance and Nd equal to 899, 
the number of WMAP temperature data points. 

The minimization of (S + eR) cannot be solved analyt- 
ically, due to the non-Gaussianity of the likelihood function, 
therefore a numerical method should be used; we found the 
iterative Newton-Raphson method suitable for our purpose. 

The m + 1-th iteration of the method is decribed by: 



5l> 



.^(S + eR) 



d.i 



(21) 



in which x™ stands for the solution i-th component relative 
to the m-th iteration. Here we have the following partial 
derivative: 



d (S + eR) 



dxi 



| £ WnN-}(\Nx m + Af) e ,(z th - z d ) e , + 
«' 

| ^(z" 1 - z d )tWtiN£QNx m + M)v{z th - z d ) e , + 

w 

i WuN-}{WJx m - d) e , + e(L*Las m )i - 

w 

\ V(Wx m - d)j(2l + l)ftW H /{\Nx m +M) 3 t - 

i Y(z th - z d ) 2 e (2e+ l)fiWu/(\Nx m +Af)t. (22) 

i 

It is intended that in the last equation z th = ln(Wa; m +Af) 
at each iteration. Furthermore the last two terms were found 
taking into account the derivative of the covariance matrix 
N _1 , assumed, just for these terms, to be diagonal and with 
diagonal elements given by: 



v -i_ (2l+l) ft 

u 2 (wx+ao?' 



(23) 



in which ft is the effective fraction of the sky computed by 
the WMAP team (Verde et al. (2003)). 

Furthermore we have approximated the Hessian matrix 
with: 



l d 2 {S + eR) 
2 dxidxj 



(W'N _1 W + eL*L) 



(24) 



where we neglected difference terms that do not substan- 
tially alter the final result and may induce numerical insta- 
bilities. At the (m + l)-th iteration step the covariance ma- 
trix is evaluated at \Nx m . Ea. l|18fl is solved for an assumed 
value of e, upon convergence the value of e is updated and 



the whole procedure is repeated until the two equations are 
both satisfied. 



2.3 Two types of errors 

With the iterative procedure we can identify the maximum 
of the a posteriori distribution of the power spectrum. Given 
the maximum and considering small deviations around the 
maximum, then the distribution function can be approxi- 
mated as Gaussian and can be fully defined by the error 
covariance matrix. To estimate the statistical significance of 
the estimated power spectrum we define two types of errors. 

According to the Bayesian method, the covariance ma- 
trix that describes the scatter of the true x around its esti- 
mator, x is given by: 



(x — x) (x — x)' 



['W t N~ 1 W + eL'Ll 



(25) 



Note that E/ provides an estimate of the deviation of the 
true x from its estimator, and it is referred to as an error of 
type I. 

Alternatively, one can ask for another measure of the 
scatter, this time of frequentist nature. Let us assume for a 
moment that the likelihood of the data is Gaussian and that 
the estimator can be expressed as: 



x = Md, 
with 

M = [W'N^W + elAl 1 W'N" 1 . 



(26) 



(27) 



Suppose that an ensemble of observers is measuring the same 
realization of the primordial power spectrum x, and that 
each observer measures and estimates it in the same way. It 
is clear that the ensemble averaged value of the estimated x 
is given by: 



v en sem b 



ble = (X ) = MVJX. 



(28) 



The variance of the scatter of the estimated x from its en- 
semble average is defined here as an error of type II, and is 
given by: 



X X enS emble ) {x Xens 



MNM 



(29) 



In other words, E// gives the scatter in the estimator of x 
around its mean value for a given realization. This corre- 
sponds to using the estimator of a; as a surrogate for the 
true model and to employing synthetic data built around it 
to study the distribution of the errors. 

Suppose that the estimation method used here, repre- 
sented by the matrix M, does not introduce any bias in x, 
then it follows that E/j provides a measure of how much x 
deviates from the true value. It will be shown that the bias 
is almost negligible for the very smooth concordance the- 
oretical spectra and that the criterion in Eg. 12011 tends to 
oversmooth the predicted signal. Consequently E/j presents 
a lower limit on the scatter of the estimated x from its the- 
oretically expected shape. 

In the next section, Monte Carlo simulations generated 
around the concordance initial spectrum are used to show 
that in the general non-Gaussian WMAP likelihood case, 
practically no bias is present and that the covariance matrix 
computed from the a posteriori distribution of the decon- 
volved samples is almost identical to E//, as expected from 
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error propagation. Furthermore the distribution very closely 
follows a Gaussian. 



2.4 Null Hypothesis Testing 

The smoothing prior necessarily induces a bias, that is 
caused by the matrix (MW) being different with respect 
to unity. In other words, the deconvolved power spectrum, 
even in absence of noise, is inevitably going to be a smoothed 
version of the true spectrum. 

However the actual concordance model power spectra 
are predicted by the simplest inflation scenarios and are 
given by very smooth shapes, namely a power law with an 
exponent n s very close to 1, that result in a very flat com- 
bination fc™ 3-1 . Therefore a smoothing prior seems reason- 
able both on theoretical grounds (since it is in line with 
inflationary predictions) and for practical purposes (since 
it offers an effective way to extract informations given the 
limitations posed by the width of the transfer function ker- 
nel). Remarkably, for extremely flat configurations the bias 
caused by the smoothing prior adopted in our method can 
be considered nearly negligible. Indeed this is what we have 
found in the Monte Carlo tests explained in the next section, 
where we considered two power law spectra with n s = 1 and 
0.97. The fact that the method has nearly negligible bias for 
configurations close to the concordance model has led us to 
consider an extremely effective Null Hypothesis Test that 
consists of the following steps: 

• the Null Hypothesis is chosen to be given, in turn, by 
the two initial spectra with n s — 1 and 0.97; 

• the CMB angular power spectra are computed from the 
initial spectra and from them random deviates are drawn to 
get synthetic data; 

• each of the Monte Carlo realizations are used as starting 
points for the reconstruction method and the deconvolved 
solutions are utilised to sample the a posteriori distribution 
of the initial spectrum around the Null Hypothesis; 

• the reconstruction from WMAP data is compared with 
the a posteriori distribution and if strong deviations are 
found, then the Null Test is considered to have failed. 

We remark that if the real primordial spectrum has de- 
viations from a smooth shape, the bias introduced by the 
method will tend to conservatively underestimate them. 

It is useful to mention that the errors of type II built 
around the Null Hypothesis spectra are expected to be 
equivalent to the errors that result from the scatter of a 
large number of Monte Carlo realizations, and indeed this is 
what we have observed. As specified in the next Section, to 
quantify the probability of how likely the Null Hypothesis 
is followed in pre-defined wavenumber ranges, we have inte- 
grated over the distribution described by the Monte Carlo 
samples. Clearly, all these considerations require that the 
important systematics have been removed from the data; 
a discussion on this delicate issue is presented in the next 
Section. 



3 RESULTS 

Since the numerical computations are rather intensive, we 
fixed the late-time cosmological parameters to the best- 



fit v alues found com b ining WMAP temperature and SDSS 
data ITegmark et alJ (12004ft : for the baryonic and matter 
fractions, Q b h 2 = 0.0231 and f2 m /i 2 = 0.149; the Hubble 
constant in units of 100 Km s" 1 Mpc -1 , h — 0.685. The opti- 
cal depth was fixed to r = 0.166, as suggested by the WMAP 
temperature and cross correlated temperature-polarization 
data, when the power spectrum is described by a power law 
with spectral index close to unity (Spergel et al. 2003). The 
range over which the transfer functions are integrated is typ- 
ically between 10 and 0.1 Mpc - , discretised in approxi- 
mately 500 bins. This choice corresponds approximatively to 
a marginalization if the precise question we are answering 
is whether there are deviations from a power law spectrum 
with spectra l index close to n 3 = 0.97, that was found in 
the best fit o flTegmark et al.l(l2004T) . while assuming a power 
spectrum defined as a power law with the spectral index as a 
free parameter. In other words, the method can be thought 
of as a consistency test for the concordance model with a 
power law power spectrum, since by fixing the cosmological 
parameters in this way, we should expect to find a power 
spectrum at least compatible with a power law of k na with 
n s — 0.97. Under these conditions the prior will tend to push 
the reconstruction towards the above power law, if it rep- 
resents the true underlying signal. While if departures from 
featureless initial conditions are found with enough statisti- 
cal significance, they may be due to a real signal or just to 
residual systematics still present in the data. We have also 
implemented the analysis with a spectral index n s — 1 to 
look for deviations from a scale free spectrum. 

As anticipated, an additional safety net is offered by 
the method: if large enough deviations truly exist, the esti- 
mate obtained with the smoothing criterion in Eg. 1201 will 
in general tend to underestimate the signal by furnishing 
an oversmoothed solution. This is a clear result of the prior 
bias that we have verified with a series of tests and is also 
deducible from a computation of the effective number of de- 
gr ees of freedom caused by s mooth ing. These are defined, as 
in ITocchini-Valentini et all ((2005) , by the trace of the ma- 
trix (WM) in analogy with ordinary regression and amount 
to about 52 for the WMAP reconstruction. Eg. 1201 1 should 
therefore be considered a rather conservative choice. 

To verify the method we performed Monte Carlo tests. 
Assuming a power law with spectral index n s — 0.97 and 
with n s = 1 as an input, we calculated the corresponding 
CMB angular power spectrum and we created random real- 
izations around it. Neglecting correlations, we sampled the 
probability distribution derived from the noise model pro- 
posed in Knox (1995): 



P(d t ) = 



y(n-2)/2 e -V/2 



Ce+Af e 2"/ 2 T(n/2) ' 
where Ci is the theoretical spectrum, n = ff(2£ + 1) and 
de + Mt, 



(30) 



V = n- 



(31) 



This distribution has the mean equal to the theoretical spec- 
trum, ( di) = Ci, and variance 



(32) 



In other words we have incorporated in the noise model the 
fit of the diagonal of the Fisher matrix N _1 obtained by 
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Figure 1. Results from 1000 Monte Carlo realizations drawn 
from a power law power spectrum, given by k" s where n a = 0.97, 
with the concordance WMAP-SDSS cosmological parameters of 
Tegmark et al.(2004). The thin smooth line in the centre is the ex- 
act initial power law and the dashed curve is the averaged recon- 
struction. Also shown are the standard deviation bands calculated 
from the scatter of the realizations around the mean. Practically 
no bias appears in the reconstruction, except at very large scales, 
where the prior starts to dominate, and however the bias is con- 
siderably smaller than the error estimation. The plot is in units 
of 2.95 X 10" 9 . 
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Figure 2. Results from 1000 Monte Carlo realizations drawn 
from a power law power spectrum, given by k" s where n s = 1, 
with the concordance WMAP-SDSS cosmological parameters of 
Tegmark et al.(2004). The thin smooth line in the centre is the ex- 
act initial power law and the dashed curve is the averaged recon- 
struction. Also shown are the standard deviation bands calculated 
from the scatter of the realizations around the mean. Practically 
no bias appears in the reconstruction, except at very large scales, 
where the prior starts to dominate, and however the bias is con- 
siderably smaller than the error estimation. The plot is in units 
of 2.95 X 10" 9 . 



the WMAP team (Verde et al. 2003) through monte carlo 
realizations of the sky that included noise, beam and cut sky 
effects. 

To speed up the computations, we considered just the 
diagonal elements of the inverse of the covariance matrix 
N , and consequently of the inverse of the matrix Q 1 in 
Eg. 11711 . and we created samples of the power separately at 
each multipole. We verified that the reconstructions oper- 
ated with the full or diagonal matrices IM" 1 and Q~ do not 
differ appreciably. Fig. and Fig. [5] show the outcome of 
the test and it is relevant to notice that the bias that results 
from the reconstruction is negligible when compared with 
type I and even type II errors. This conclusion supports our 
consideration of the reliability of the less conservative type 
II errors, when utilised to detect the departure from the 
flat and concordance power law models. We also noted that 
type II errors are almost identical to the standard deviation 
uncertainties computed from the Monte Carlo samples; this 
provides a useful check for the method. 

For all the reconstructions from the Monte Carlo re- 
alizations, we fix the parameter e to the value that results 
from deconvolving the WMAP data according to Eo. l20ll . 
We have also checked that when choosing e following Eg. 1201 1 
for each individual Monte Carlo sample built from an ini- 
tial spectrum with features added, the resulting mean recon- 
struction was a conservative smoothed version of the starting 
one. 

The analysis of the Monte Carlo realizations validates 
that the algorithm used here indeed produces an almost un- 
biased estimator of power spectrum, when very flat spectra 
are used as input. If the power spectrum has real features, 
then some bias will be inevitably present in the estimator, 
but it will act in order to conservatively underestimate the 
features. 

We verified that the Monte Carlo reconstructions on 
all scales are well described by a Gaussian distribution that 
follows type II errors. 

In Fig. [3] our most important result is plotted, namely 



the estimated power spectrum from the WMAP temper- 
ature data (Bennett et al. 2003) at large scales together 
with the two type of errors. In Fig. 2] we show the WMAP 
reconstruction compared with the Monte Carlo results de- 
scribed previously for the power law spectrum case. The 
figure can provide intuition of how unlikely in some wave 
number ranges the reconstruction could resemble a statisti- 
cal random deviation from the fiducial spectrum, and depicts 
a graphical version of the Null Hypothesis Test described in 
the previous Section. 

Assuming the flat spectrum to define the fiducial initial 
conditions, we find evidence for the three following potential 
deviations from a power spectrum without features consid- 
ering type I (type II) errors: a cutoff for the largest scales 
0.0001 < k < 0.001 Mpc" 1 at 1 (2) sigma c.L; a dip at 
0.001 < k < 0.003 Mpc" 1 at 1.5 (3) sigma c.L; a bump a 
0.003 < k < 0.005 Mpc" 1 at 1 (2) sigma c.L. While with 
the power law as the fiducial spectrum, we find, considering 
type I (type II) errors: a cutoff at 0.0001 < k < 0.001 Mpc" 1 
at 1.5 (2.5) sigma c.L; a dip at 0.001 < k < 0.003 Mpc" 1 at 
1.5 (3) sigma c.L; a bump a 0.003 < k < 0.005 Mpc" 1 at 1 
(2) sigma c.L. 

We now provide another way to estimate the deviations, 
that is practically equivalent to integrating over the distri- 
bution defined by the errors of type II built around the Null 
Hypothesi s spectra in defined ranges in wavenumber space. 
Following iKoeo et"aT1 l)2004afl , we define a measure of the 
distance between a generic spectrum (P 9 ) and the fiducial 
one (P ) in the wavenumber range limited by k\ and k^\ 



D(k u k 2 ) 



fci 



k2 dk (P g (k) P f (k) 



P 9 {h) _ P f (kj) 



(33) 



We calculated the distance of all the Monte Carlo samples 
and the WMAP reconstruction with respect to the fiducial 
spectrum and we counted how many realizations had a dis- 
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Figure 3. Reconstructed initial power spectrum in units of 
2.95 X 10 — 9 from WMAP data with the concordance model cos- 
mological parameters. The irregular middle curve represent the 
mean and the surrounding curves that delimit the light shaded 
band are obtained by summing and subtracting the square root of 
the diagonal elements of the covariance matrix for type II errors. 
The borders of the larger dark shaded region are given by sum- 
ming and subtracting the square root of the diagonal elements of 
the covariance matrix for type I errors. The smooth line is a ref- 
erence power law power spectrum given by k" s , where n s = 0.97. 
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Figure 4. Results from 1000 Monte Carlo realizations drawn 
from a power law power spectrum, given by k" s where n s = 0.97, 
with the concordance WMAP-SDSS cosmological parameters of 
Tegmark et al.(2004). The thin smooth line in the centre is the ex- 
act initial power law and the dashed curve is the averaged recon- 
struction. Also shown are the standard deviation bands calculated 
from the scatter of the realizations around the mean. The thick 
continuous line, the reconstruction from WMAP data, is also plot- 
ted to show how likely it is that the deconvolution could represent 
a random deviate of the power law fiducial model, assumed to be 
our Null Hypothesis. The plot is in units of 2.95 X 10~ 9 . 



tance larger than the one computed for the WMAP data de- 
convolution. In this way we could quantify the probabilities 
of exceeding the observed deviations in WMAP data respect 
to the fiducial models. We detected the following potential 
departures from a scale free (spectral index n s = 0.97) ini- 
tial spectrum: a cut-off at 0.0001 < k < 0.001 Mpc" 1 at 
79.5% (92%), a dip at 0.001 < k < 0.003 Mpc" 1 at 87.2% 
(98%) and a bump at 0.003 < k < 0.004 Mpc -1 at 90.3% 
(55.5%). 

We add that the pri or, similarly to the Weiner filter (see 
e.g. lZaroubi et all l)l995h '). helps in retrieving useful informa- 
tion at the crossover between the signal and noise dominated 
regimes, if combined with Monte Carlo testing. Furthermore 
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Figure 5. The thick dashed smooth line represents the re- 
constructed power spectrum from WMAP data convolved back 
to multipole space, while the power law with spectral index 
n s = 0.97, best fit of the ACDM cosmology, is plotted as a con- 
tinuous thick curve. The WMAP data are also shown as binned 
data points, with the error bars computed assuming the power 
law spectrum. 



methods such as ours offer the possibility to penetrate the 
informational barrier built up by cosmic variance. 

We plot in Fig. the WMAP data together with the 
convolved estimated power spectrum in multipole space to 
illustrate the effect of smoothing. That interesting features 
might have been present in the power spectrum could have 
already been guessed from an inspection of the WMAP data 
confronted with a power law power spectrum ACDM pre- 
diction. However how the information (including the exper- 
imental uncertainties) in multipole space was going to be 
transferred in detail to wave number space could not have 
been predicted reliably without an analysis that had taken 
properly into account the complicated transfer functions. 
Furthermore, the processed information is now in a format 
that permits a direct comparison with the theoretical pre- 
dictions from inflation. 

Even though we compute the full covariance matrices, 
in the plots we show the errors given by the square root of 
the inverse diagonal elements of the relative covariance ma- 
trix for illustrational purposes. However it is important to 
assess if the deviations detected and their magnitude take 
their origin from error correlation. We implemented some 
tests to check this relevant issue. In turn, we replaced the 
WMAP multipoles that correspond to the scales of the three 
detected features with the concordance model predictions 
without any noise. From the reconstructions in Fig. [fj] wc 
deduce that the evidence for the features is not the product 
of correlation, since their significance is only slightly wors- 
ened. It is interesting to note that, as expected, the low 
values of the dipole and the quadrupole are the principal 
cause for the observed cutoff. Furthermore the wiggles in 
the WMAP data around the multipoles of 20 and 40 show 
up clearly as the bump and the dip in our reconstruction. 

We show Fig. |7| how the reconst r uction changes if we 
adopt the analysis of iTeemark et all (l2005tl . based on an 
all-sky foreground cleaned map derived from WMAP data. 
For representative purposes, we utilise the WMAP data like- 
lihood with the mean values of the quadrupole and octopole 
lifted to respectively 202 and 870 /iK 2 , rather than using the 
original WMAP values of 123 and 612 fiK 2 . It can be seen 
that, with respect to the concordance power law spectrum 
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Figure 6. The figure shows the effects of correlation and sheds 
light on the multipole space origin of the features detected. To 
serve such a purpose, we replaced some WMAP data points with 
the values of prediction of the ACDM concordance model with a 
power law power spectrum given by k" s , where n s = 0.97. The 
dotted curve is the WMAP reconstruction and the error bands 
are derived from type II errors (see text). Going from top left 
and clockwise the power was modified at the following multipoles: 
i = 2; I = 2, 3; 15 < I < 26: 26 < I < 56. 

with spectral index n 3 — 0.97, a cut-off detection is still 
present at about 2a in frequentist terms, while the Bayesian 
estimate finds the x of the Null Hypothesis lcr away from the 
estimate x. We report that some other of the post- WMAP 
papers on the estimation of the quadrupole point towards 
a somewhat larger value respect to the original WMAP one 
(e.g. Efstathiou 2004 and Slosar et al. 2004). For example 
Slosar et al. (2004) find a broad quadrupole likelihood func- 
tion with a tail that favours greater values, as a result of 
having marginalized over the foreground templates. This dif- 
fers from the WMAP analysis that consisted in fitting and 
subtracting out the templates from the data. Clearly both 
the estimated mean and errors of the reconstructed initial 
power spectrum would be affected by the modification of the 
CMB dataset. 

We also present in Fig.|S]a reconstruction with the best 
fit parameters found in Spergel et al. (2003) , when consider- 
ing a power law spectrum: Qbh 2 = 0.023 and f2 m /i 2 = 0.13; 
the Hubble constant in units of 100 Kms" 1 Mpc" 1 , h — 0.68; 
the optical depth, r = 0.1. The main effect with respect to 
the model considered previously is a reduction of the cut-off 
significance due to a smaller optical depth, as the main ef- 
fect would be to reduce the power at intermediate and small 
scales while leaving the power unaltered at large scales. A 
conspiracy between a low optical depth and enhanced power 
in the quadrupole and octopole could in principle consider- 
ably decrease the detection of a cut-off; such measurements 
should therefore be followed closely in future experiments. 



4 DISCUSSION 

Given the temperature data from the WMAP experiment we 
have reconstructed at high resolution the primordial power 
spectrum at large scales by adopting an iterative smoothing 
algorithm, built to deal with the non-gaussianity of the data 
errors caused mainly by cosmic variance. We have fixed the 
paramet ers to the best fi t late-t ime cosmological parameters 
found bv lTeemark et all ll2004l) in order to search for devia- 
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Figure 7. The result of considering quadrupole and octopole 
mean values raised to the Tegmark et al. (2003) all-sky foreground 
subtracted map determination. The smooth curve represents a 
power-law power spectrum given by fc Us , where n s = 0.97, while 
the dotted line is the deconvolution from the original WMAP 
data. The continuous curve is the reconstruction with the new 
data, the error bands come from the type I (darker) and type II 
(lighter) errors. The smoothing parameter e was fixed according 
to Ea. 1201 using the new data set. 
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Figure 8. We show the reconstruction for the cosmological late- 
time parameters given in Spergel et al. (2003), where an optical 
depth t = 0.1 is considered. The smooth curve represents a power 
law power spectrum given by fc™ s , where n s = 0.97, while the 
dotted line is the deconvolution from the original WMAP data. 
The continuous curve is the reconstruction with the new parame- 
ters, the error bands coming from the type I (darker) and type II 
(lighter) errors. The smoothing parameter e was fixed according 
to Eg. 1201 . taking into account the new transfer function. The 
smaller optical depth causes a reduction in the statistical signifi- 
cance of the cutoff. 

tions from the concordance featurless power law power spec- 
trum, with spectral index n s approximately equal to 0.97. 
We also implemented the analysis with a scale free spec- 
trum. The horizon scale for the model we are considering 
corresponds to the wavenumber kh = 4.52 x 10~ 4 Mpc" 1 . By 
creating and reconstructing from Monte Carlo realizations 
of CMB data, we have verified that our proposed estima- 
tor is effectively unbiased if a scale free or power law with 
spectral index n s — 0.97 power spectra are assumed. This 
has led us to the detection of the following three deviations 
from a scale free (spectral index n s — 0.97) initial spectrum: 
a cut-off at 0.0001 < k < 0.001 Mpc" 1 at 79.5% (92%), a dip 
at 0.001 < k < 0.003 Mpc" 1 at 87.2% (98%) and a bump at 
0.003 < k < 0.004 Mpc" 1 at 90.3% (55.5%). 
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We have also derived using two alternative ways a mea- 
sure of the local uncertainties in the estimates of deviations 
and effective location of the initial power spectrum, labeled 
by errors of type I and II. The frequentist analysis finds the 
low k cutoff of the estimated power spectrum to be at about 
2.5o/i away from the n s = 0.97 model, while in the Bayesian 
analysis the model is about 1.5<r/ away from the estimated 
spectrum. 

The convergence in the observed low power at large 
scales by the WMAP and COBE DMR experiments is cer- 
tainly interesting and merits attention. Since further high 
quality large-scale temperature and polarization CMB data 
are expected in the near future, for example from the next 
data releases of the WMAP team and the future Planck mis- 
sion, any method that could shed light on the nature of the 
primordial perturbations will prove to be extremely useful 
for improving our understanding of the mechanism that is 
responsible for the generation of cosmic structure. 
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